Reproducible Analysis
Fully reproducible analysis stream
Lombardo et al., eLife
# Load necessary libraries
library(easypackages)
libraries("nlme","matlabr","ggplot2","multcomp","readxl",
"heplots","effsize","reshape2","plyr","psych",
"MASS","here","rsample")
# paths
codedir = here("code")
tidydatadir = here("data","tidy")
resultdir = here("results")
plotdir = here("plots")
# allows for selecting n colors from color wheel like ggplot2
source(file.path(codedir,"get_ggColorHue.R"))
# for running clinical trajectory analyses
source(file.path(codedir,"functions4trajanalysis.R"))
# Set options and other stuff
options(stringsAsFactors = FALSE)
options(tibble.print_max = Inf)
options(matlab.path = "/Applications/MATLAB_R2018b.app/bin")
fontSize = 10
dotSize = 3
fdr_thresh = 0.05
RUNMATLAB = FALSERead in pheno data
Dfmri = read_excel(file.path(tidydatadir,"final_allETrsfMRIsubs_phenodata04_ASDTDLDDDSIB.xlsx"))
Dfmri$subjectId = factor(Dfmri$subjectId)
Dfmri$subgrp2 = factor(Dfmri$subgrp2)
Dfmri$subgrp2 = factor(Dfmri$subgrp2,levels(Dfmri$subgrp2)[c(2,4,3,6,5,1)])
colours2use = get_ggColorHue(6)Analyze demographics
demoVar_sub = Dfmri[,c("subgrp2","scan_age","sex",
"ET1_Age","meanFD",
"meanDVARSraw","meanDVARSwavelet")]
describeBy(demoVar_sub, group = "subgrp2")##
## Descriptive statistics by group
## subgrp2: GeoASD
## vars n mean sd median trimmed mad min max range
## subgrp2* 1 16 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00
## scan_age 2 16 29.93 8.72 28.04 30.06 11.50 14.16 43.79 29.63
## sex* 3 16 NaN NA NA NaN NA Inf -Inf -Inf
## ET1_Age 4 16 28.38 7.78 27.00 28.43 7.41 15.00 41.00 26.00
## meanFD 5 16 0.06 0.02 0.06 0.06 0.02 0.03 0.11 0.08
## meanDVARSraw 6 16 6.87 1.14 6.99 6.94 1.26 4.45 8.38 3.93
## meanDVARSwavelet 7 16 5.31 0.87 5.33 5.32 0.79 3.79 6.72 2.93
## skew kurtosis se
## subgrp2* NaN NaN 0.00
## scan_age 0.03 -1.20 2.18
## sex* NA NA NA
## ET1_Age 0.07 -1.12 1.94
## meanFD 0.60 -1.14 0.01
## meanDVARSraw -0.55 -0.85 0.28
## meanDVARSwavelet -0.04 -1.15 0.22
## --------------------------------------------------------
## subgrp2: nonGeoASD
## vars n mean sd median trimmed mad min max range
## subgrp2* 1 62 2.00 0.00 2.00 2.00 0.00 2.00 2.00 0.00
## scan_age 2 62 29.37 8.35 30.47 29.65 9.04 12.35 44.06 31.70
## sex* 3 62 NaN NA NA NaN NA Inf -Inf -Inf
## ET1_Age 4 62 26.81 8.35 26.00 26.60 7.41 12.00 44.00 32.00
## meanFD 5 62 0.10 0.12 0.06 0.07 0.02 0.03 0.93 0.90
## meanDVARSraw 6 62 7.00 1.61 6.70 6.86 1.00 4.14 15.17 11.03
## meanDVARSwavelet 7 62 5.34 0.79 5.27 5.31 0.75 3.65 7.30 3.65
## skew kurtosis se
## subgrp2* NaN NaN 0.00
## scan_age -0.35 -0.84 1.06
## sex* NA NA NA
## ET1_Age 0.24 -0.74 1.06
## meanFD 5.32 32.56 0.02
## meanDVARSraw 2.11 8.55 0.20
## meanDVARSwavelet 0.32 0.00 0.10
## --------------------------------------------------------
## subgrp2: LD_DD
## vars n mean sd median trimmed mad min max range
## subgrp2* 1 15 3.00 0.00 3.00 3.00 0.00 3.00 3.00 0.00
## scan_age 2 15 25.12 7.97 23.95 24.90 9.35 13.37 39.75 26.38
## sex* 3 15 NaN NA NA NaN NA Inf -Inf -Inf
## ET1_Age 4 11 19.36 4.15 20.00 19.44 4.45 13.00 25.00 12.00
## meanFD 5 15 0.10 0.05 0.08 0.09 0.04 0.03 0.19 0.15
## meanDVARSraw 6 15 7.41 2.14 6.82 7.08 1.10 5.16 14.00 8.84
## meanDVARSwavelet 7 15 5.55 1.22 5.39 5.37 0.78 4.11 9.33 5.21
## skew kurtosis se
## subgrp2* NaN NaN 0.00
## scan_age 0.28 -1.06 2.06
## sex* NA NA NA
## ET1_Age -0.26 -1.44 1.25
## meanFD 0.60 -1.17 0.01
## meanDVARSraw 1.88 3.22 0.55
## meanDVARSwavelet 1.80 3.28 0.32
## --------------------------------------------------------
## subgrp2: TypSibASD
## vars n mean sd median trimmed mad min max range
## subgrp2* 1 16 4.00 0.00 4.00 4.00 0.00 4.00 4.00 0.00
## scan_age 2 16 26.74 9.38 27.89 26.51 9.74 12.52 44.09 31.57
## sex* 3 16 NaN NA NA NaN NA Inf -Inf -Inf
## ET1_Age 4 14 19.79 6.20 19.50 19.42 8.15 13.00 31.00 18.00
## meanFD 5 16 0.08 0.04 0.07 0.08 0.03 0.03 0.18 0.15
## meanDVARSraw 6 16 6.82 1.57 6.47 6.69 0.82 4.64 10.77 6.13
## meanDVARSwavelet 7 16 5.28 0.88 5.14 5.22 0.49 3.91 7.44 3.53
## skew kurtosis se
## subgrp2* NaN NaN 0.00
## scan_age 0.02 -1.01 2.35
## sex* NA NA NA
## ET1_Age 0.41 -1.31 1.66
## meanFD 0.75 -0.63 0.01
## meanDVARSraw 0.99 0.31 0.39
## meanDVARSwavelet 0.90 0.39 0.22
## --------------------------------------------------------
## subgrp2: TD
## vars n mean sd median trimmed mad min max
## subgrp2* 1 55 5.00 0.00 5.00 5.00 0.00 5.00 5.00
## scan_age 2 55 29.61 10.14 30.92 29.70 11.69 13.17 47.93
## sex* 3 55 NaN NA NA NaN NA Inf -Inf
## ET1_Age 4 55 23.07 9.07 22.00 22.29 11.86 12.00 45.00
## meanFD 5 55 0.11 0.11 0.07 0.08 0.03 0.04 0.59
## meanDVARSraw 6 55 7.00 1.44 6.85 6.87 1.52 4.63 11.32
## meanDVARSwavelet 7 55 5.19 0.61 5.11 5.17 0.66 3.89 6.58
## range skew kurtosis se
## subgrp2* 0.00 NaN NaN 0.00
## scan_age 34.76 -0.21 -1.12 1.37
## sex* -Inf NA NA NA
## ET1_Age 33.00 0.57 -0.68 1.22
## meanFD 0.55 2.66 7.20 0.01
## meanDVARSraw 6.69 0.79 0.34 0.19
## meanDVARSwavelet 2.69 0.25 -0.48 0.08
## --------------------------------------------------------
## subgrp2: ASDnoET
## vars n mean sd median trimmed mad min max range
## subgrp2* 1 31 6.00 0.00 6.00 6.00 0.00 6.00 6.00 0.00
## scan_age 2 31 29.69 8.88 30.03 29.81 11.50 13.21 43.63 30.42
## sex* 3 31 NaN NA NA NaN NA Inf -Inf -Inf
## ET1_Age 4 0 NaN NA NA NaN NA Inf -Inf -Inf
## meanFD 5 31 0.09 0.08 0.06 0.07 0.02 0.03 0.41 0.37
## meanDVARSraw 6 31 6.98 1.46 6.91 6.93 1.05 4.00 11.63 7.63
## meanDVARSwavelet 7 31 5.24 0.76 5.11 5.26 0.71 3.57 6.56 2.99
## skew kurtosis se
## subgrp2* NaN NaN 0.00
## scan_age -0.05 -1.22 1.59
## sex* NA NA NA
## ET1_Age NA NA NA
## meanFD 2.42 5.45 0.01
## meanDVARSraw 0.73 1.74 0.26
## meanDVARSwavelet -0.11 -0.58 0.14Scan Age ANOVA
mod2use = lm(scan_age ~ subgrp2, data = demoVar_sub)
anova(mod2use)## Analysis of Variance Table
##
## Response: scan_age
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp2 5 366 73.208 0.8914 0.4879
## Residuals 189 15522 82.127Eye tracking age ANOVA
mod2use = lm(ET1_Age ~ subgrp2, data = demoVar_sub)
anova(mod2use)## Analysis of Variance Table
##
## Response: ET1_Age
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp2 4 1283.4 320.84 4.7742 0.001174 **
## Residuals 153 10282.0 67.20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sex by Group Table for all individuals
subgrp_sex_tab = table(demoVar_sub$subgrp2, demoVar_sub$sex, exclude = "NA")
subgrp_sex_chisq_res = chisq.test(subgrp_sex_tab)
subgrp_sex_tab##
## F M
## GeoASD 5 11
## nonGeoASD 13 49
## LD_DD 5 10
## TypSibASD 8 8
## TD 18 37
## ASDnoET 4 27subgrp_sex_chisq_res##
## Pearson's Chi-squared test
##
## data: subgrp_sex_tab
## X-squared = 9.8871, df = 5, p-value = 0.0785Analyze Head Motion Measures for Group Differences and make plots
Will look at mean framewise displacement (mean FD in mm), as well as DVARS measurements before and after wavelet denoising, just to show the impact wavelet denoising has on removing substantial amounts of artifact from the data.
# Mean Framewise Displacement Plot
p = ggplot(data = Dfmri, aes(x = subgrp2, y = meanFD, colour = subgrp2))
p = p + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA)
p = p + scale_colour_manual(values = colours2use) + ylab("Mean Framewise Displacement (mm)")
ggsave(filename = file.path(plotdir,"meanFDplot.pdf"))
p
Make Mean DVARS Plot that shows DVARS before and after wavelet denoising
Dsub = melt(as.data.frame(Dfmri), id = c("subjectId","subgrp2"),
measure = c("meanDVARSraw","meanDVARSwavelet"))
Dsub$value = as.numeric(Dsub$value)
Dsub$variable = factor(Dsub$variable)
Dsub$variable = revalue(Dsub$variable, c("meanDVARSraw"="Before Wavelet Denoising",
"meanDVARSwavelet"="After Wavelet Denoising"))
p = ggplot(data = Dsub, aes(x = subgrp2, y = value, colour = subgrp2)) + facet_grid(. ~ variable)
p = p + geom_jitter(size = dotSize)+
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA)
p = p + scale_colour_manual(values = colours2use) + ylab("Mean DVARS (%)")
ggsave(filename = file.path(plotdir,"meanDVARSplot.pdf"))
p
Run ANOVAs on measures of head motion (meanFD)
# Mean FD
lm_formula = as.formula(sprintf("%s ~ %s","meanFD","subgrp2"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dfmri, na.action = na.omit)))
anova(mod2use)## Analysis of Variance Table
##
## Response: meanFD
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp2 5 0.03298 0.0065955 0.6759 0.6422
## Residuals 189 1.84420 0.0097577Mean DVARS before and after wavelet denoising
# Before
lm_formula = as.formula(sprintf("%s ~ %s","meanDVARSraw","subgrp2"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dfmri, na.action = na.omit)))
anova(mod2use)## Analysis of Variance Table
##
## Response: meanDVARSraw
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp2 5 3.38 0.67586 0.2811 0.9231
## Residuals 189 454.45 2.40450# After
lm_formula = as.formula(sprintf("%s ~ %s","meanDVARSwavelet","subgrp2"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dfmri, na.action = na.omit)))
anova(mod2use)## Analysis of Variance Table
##
## Response: meanDVARSwavelet
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp2 5 1.808 0.36151 0.5716 0.7217
## Residuals 189 119.530 0.63244Make plot of eye tracking data
This plot is from eye tracking data from just the subjects who also have rsfMRI data available.
Dsub = subset(Dfmri, !Dfmri$subgrp2=="ASDnoET", select = 1:ncol(Dfmri))
Dsub$subgrp2 = factor(Dsub$subgrp2)
Dsub$Dx = factor(Dsub$Dx)
Dsub$Dx = factor(Dsub$Dx,levels(Dsub$Dx)[c(1,2,4,3)])
c2use = c(colours2use[1], colours2use[5], colours2use[2:4])
xLabel = "Group"
yLabel = "Fixation Geometric (%)"
p = ggplot(data = Dsub, aes(x = Dx, y = Percent_Fixation_Geometric/100, colour = subgrp2))
p = p + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
guides(colour = FALSE)
p = p + scale_y_continuous(limits = c(0,1),breaks = round(seq(from = 0, to = 1, by = 0.1),digits=2))
p = p + geom_hline(yintercept = 0.69, linetype = 2) +
scale_colour_manual(values = c2use) +
xlab(xLabel) + ylab(yLabel) +
theme(axis.text.x = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+10,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+10,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+10,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+10,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,"eyetracking_geoFix_plot.pdf"))
p
Read in data for longitudinal clinical trajectory analyses
# Read in longitudinal Mullen, Vineland, and ADOS summary score data
Dlw = read_excel(file.path(tidydatadir,"LW_Report_12012016.xlsx"), na = "NULL")Pull out data
# masks to pull out groups
na_mask = is.na(Dlw[,1])
asd_mask = is.element(Dlw$dxCode, ASDlabels) & !na_mask
td_mask = is.element(Dlw$dxCode, TDlabels) & !na_mask
dd_mask = is.element(Dlw$dxCode, DDlabels) & !na_mask
ld_mask = is.element(Dlw$dxCode, LDlabels) & !na_mask
other_mask = is.element(Dlw$dxCode, OTHERlabels) & !na_mask
asdfeat_mask = is.element(Dlw$dxCode, ASDFEATlabels) & !na_mask
typsib_mask = is.element(Dlw$dxCode, TYPSIBlabels) & !na_mask
# column labels of interesting data to extract
allLabels = c("VageMo","ComTotal_DomStd","DlyTotal_DomStd","SocTotal_DomStd",
"MtrTotal_DomStd","AdapBehav_DomStd",
"AageMo","CoSoTot","RRTot","CoSoTotRRTot",
"MageMo","VRT","VR_Raw","FMT","FM_Raw",
"RLT","RL_Raw","ELT","EL_Raw")
# grab the data for each group
# ASD
asd_df = grabGroupData(D = Dlw, submask = asd_mask, colLabels = allLabels, grpLabel = "ASD")
Det = read.delim(file.path(tidydatadir,"LW_Report_for_ET_N937.txt"),na.strings = c("NA","NULL"))
asd_df = getETsubgrp2(tidy_df = asd_df, full_df = Det, Dx = "ASD")
asd_df = subset(asd_df, asd_df$ETsubgrpDx!="ASD", select = 1:ncol(asd_df))
asd_df = cleanAgeErrors(asd_df)
asd_df$subjectId = factor(asd_df$subjectId)
nongeo_df = subset(asd_df, is.element(asd_df$subjectId,Dfmri$subjectId) & asd_df$ETsubgrpDx=="nonGeo ASD",
select = 1:ncol(asd_df))
geo_df = subset(asd_df, asd_df$ETsubgrpDx=="Geo ASD",
select = 1:ncol(asd_df))
fmri_df = rbind(nongeo_df, geo_df)
fmri_df$ETsubgrpDx = factor(fmri_df$ETsubgrpDx)
fmri_df$p2f2 = factor(fmri_df$p2f2)
fmri_df$Dx = factor(fmri_df$Dx)
# set general stuff for the plot
cols2use = get_ggColorHue(6)
cols2use = c(cols2use[1], cols2use[5])
# what data and what subgroup do you want to analyze
df2use = fmri_df
df2use$ETsubgrpDx = as.factor(df2use$ETsubgrpDx)
subgrp_var = "ETsubgrpDx"Mullen Receptive Language
modelType = "linear"
xLabel = "Age (Months)"
plot_dots = TRUE
plot_lines = TRUE
ci_band = TRUE
dot_alpha = 5/10
line_alpha = 5/10
band_alpha = 4/10
standardize = TRUE
x_var = "MageMo"
xLimits = c(5,55)
yLimits = c(0,75)
yLabel = "Mullen Receptive Language T-Score (RL)"
y_var = "RLT"
RLT_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(RLT_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 165 0.026293 0.8714
## MageMo 1 165 0.714530 0.3992
## ETsubgrpDx 1 117 21.927787 <.0001
## MageMo:ETsubgrpDx 1 165 5.223953 0.0236# change the coloring of the groups to match other figures
p = RLT_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Mullen_%s_traj_plot.pdf",y_var)))
p
Mullen Expressive Language
yLabel = "Mullen Expressive Language T-Score (EL)"
y_var = "ELT"
ELT_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(ELT_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 165 0.243046 0.6227
## MageMo 1 165 3.417672 0.0663
## ETsubgrpDx 1 117 21.185912 <.0001
## MageMo:ETsubgrpDx 1 165 1.074913 0.3014# change the coloring of the groups to match other figures
p = ELT_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Mullen_%s_traj_plot.pdf",y_var)))
p
Mullen Visual Reception
yLabel = "Mullen Visual Reception T-Score (VR)"
y_var = "VRT"
VRT_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VRT_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 165 1.297810 0.2563
## MageMo 1 165 20.157709 <.0001
## ETsubgrpDx 1 117 7.910200 0.0058
## MageMo:ETsubgrpDx 1 165 5.639771 0.0187# change the coloring of the groups to match other figures
p = VRT_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Mullen_%s_traj_plot.pdf",y_var)))
p
Mullen Fine Motor
yLabel = "Mullen Fine Motor T-Score (FM)"
y_var = "FMT"
FMT_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(FMT_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 165 0.14507 0.7038
## MageMo 1 165 76.64519 <.0001
## ETsubgrpDx 1 117 15.33978 0.0002
## MageMo:ETsubgrpDx 1 165 0.90434 0.3430# change the coloring of the groups to match other figures
p = FMT_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Mullen_%s_traj_plot.pdf",y_var)))
p
Vineland Communication
x_var = "VageMo"
xLimits = c(0,55)
yLimits = c(25,125)
yLabel = "Vineland Communication Standard Score"
y_var = "ComTotal_DomStd"
VineComm_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VineComm_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 185 3.439942 0.0652
## VageMo 1 185 1.514647 0.2200
## ETsubgrpDx 1 120 17.142872 0.0001
## VageMo:ETsubgrpDx 1 185 0.179004 0.6727# change the coloring of the groups to match other figures
p = VineComm_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Vineland_%s_traj_plot.pdf",y_var)))
p
Vineland Socialization
yLabel = "Vineland Socialization Standard Score"
y_var = "SocTotal_DomStd"
VineSoc_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VineSoc_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 185 0.00000 0.9985
## VageMo 1 185 35.98014 <.0001
## ETsubgrpDx 1 120 9.91304 0.0021
## VageMo:ETsubgrpDx 1 185 3.23199 0.0738# change the coloring of the groups to match other figures
p = VineSoc_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Vineland_%s_traj_plot.pdf",y_var)))
p
Vineland Daily Living
yLabel = "Vineland Daily Living Standard Score"
y_var = "DlyTotal_DomStd"
VineDly_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VineDly_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 185 0.133603 0.7151
## VageMo 1 185 15.883371 0.0001
## ETsubgrpDx 1 120 12.239522 0.0007
## VageMo:ETsubgrpDx 1 185 3.699051 0.0560# change the coloring of the groups to match other figures
p = VineDly_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Vineland_%s_traj_plot.pdf",y_var)))
p
Vineland Motor
yLabel = "Vineland Motor Standard Score"
y_var = "MtrTotal_DomStd"
VineMtr_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VineMtr_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 185 0.62540 0.4301
## VageMo 1 185 46.10245 <.0001
## ETsubgrpDx 1 120 4.89719 0.0288
## VageMo:ETsubgrpDx 1 185 1.60509 0.2068# change the coloring of the groups to match other figures
p = VineMtr_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Vineland_%s_traj_plot.pdf",y_var)))
p
Vineland Adaptive Behavior
yLabel = "Vineland Adaptive Behavior Standard Score"
y_var = "AdapBehav_DomStd"
VineAdapBehav_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(VineAdapBehav_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 185 0.229060 0.6328
## VageMo 1 185 13.748281 0.0003
## ETsubgrpDx 1 120 13.889196 0.0003
## VageMo:ETsubgrpDx 1 185 1.967818 0.1624# change the coloring of the groups to match other figures
p = VineAdapBehav_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("Vineland_%s_traj_plot.pdf",y_var)))
p
ADOS Social Affect
x_var = "AageMo"
xLimits = c(0,55)
yLimits = c(0,25)
yLabel = "ADOS Social Affect"
y_var = "CoSoTot"
ADOSCoSo_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(ADOSCoSo_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 181 0.975626 0.3246
## AageMo 1 181 0.549231 0.4596
## ETsubgrpDx 1 120 9.074773 0.0032
## AageMo:ETsubgrpDx 1 181 3.418836 0.0661# change the coloring of the groups to match other figures
p = ADOSCoSo_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("ADOS_%s_traj_plot.pdf",y_var)))
p
ADOS Repetitive Restricted Behavior
yLabel = "ADOS Repetitive Restricted Behavior"
y_var = "RRTot"
yLimits = c(0,15)
ADOSRRB_all_quad = spaghettiPlot(df = df2use, x_var = x_var, y_var = y_var,
subgrp_var = subgrp_var, xLabel = xLabel, modelType = modelType,
yLabel = yLabel, fname2save = NULL,
plot_dots = plot_dots, plot_lines = plot_lines, ci_band = ci_band,
dot_alpha = dot_alpha, line_alpha = line_alpha, band_alpha = band_alpha,
xLimits = xLimits, yLimits = yLimits, standardize = standardize)
anova(ADOSRRB_all_quad$lme_model)## numDF denDF F-value p-value
## (Intercept) 1 181 0.211298 0.6463
## AageMo 1 181 0.131105 0.7177
## ETsubgrpDx 1 120 5.038428 0.0266
## AageMo:ETsubgrpDx 1 181 1.031324 0.3112# change the coloring of the groups to match other figures
p = ADOSRRB_all_quad$p
p = p + scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use) +
theme(axis.text.x = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("ADOS_%s_traj_plot.pdf",y_var)))
p
Compute partial correlations in MATLAB and set up data to run tests
Use R to run MATLAB code to estimate partial correlations (Tikhonov regularisation) with FSLNets code. Relies on MATLAB function estimateConnectivity.m being placed in your code directory.
if (RUNMATLAB){
code2run = sprintf("cd %s; estimateConnectivity;",codedir)
res = run_matlab_code(code2run)
}
dfname = file.path(tidydatadir,"partialCorDataASDTDLDDDSIB_ridge_lambda1.txt")
D = read.delim(dfname, na.strings = c("NA", "NaN", " "))Split up data into subsets for further analyses.
# make sure variables are factors
D$subjectId = factor(D$subjectId)
D$sex = factor(D$sex)
D$subgrp = factor(D$subgrp)
D$subgrp = factor(D$subgrp,levels(D$subgrp)[c(2,4,3,6,5,1)])
D$CaseControl = factor(D$CaseControl)
# find variables names for the IC-pairs
vars2use = colnames(D)[8:ncol(D)]
D$CaseControl2 = as.character(D$CaseControl)
D$CaseControl2[D$subgrp=="LD_DD"] = "LD_DD"
D$CaseControl2[D$subgrp=="TypSibASD"] = "TypSibASD"
D$CaseControl2 = factor(D$CaseControl2)
# grab subsets of data for expt 2
Dexp2_all = subset(D, is.element(D$subgrp,c("GeoASD","nonGeoASD","LD_DD","TypSibASD","TD")),
select = 1:ncol(D))
Dexp2_all$subgrp = factor(Dexp2_all$subgrp)
colours2use = get_ggColorHue(6)
colours2use = c(colours2use[1],colours2use[5],colours2use[2:4])
Dtmp = merge(Dexp2_all,Dfmri, by = "subjectId")
Dtmp$subgrp = factor(Dtmp$subgrp.x)
Dtmp$sex = factor(Dtmp$sex.x)
Dtmp$scan_age = Dtmp$scan_age.xRun ANOVAs on each IC-pair on CaseControl, Sex, Scan Age model
colnames2use = c("df1_all","df2_all",
"Fstat_all","pval_all","etasq_all","fdr_all",
"df1_subtype","df2_subtype",
"fstat_subtype","pval_subtype","etasq_subtype","fdr_subtype",
"fstat_subtypeDVARScov","pval_subtypeDVARScov","fdr_subtypeDVARScov",
"AIC_nosubtype","AIC_subtype","AIC_delta")
aovres = data.frame(matrix(nrow = length(vars2use),ncol =length(colnames2use)))
colnames(aovres) = colnames2use
rownames(aovres) = vars2use
for (i in 1:length(vars2use)) {
y_var = vars2use[i]
# construct linear model for ASD vs TD vs LD/DD vs TD ASDSib
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"CaseControl2","sex","scan_age"))
mod2use = eval(substitute(lm(formula = lm_formula, data = D, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
# extract F-stat and pvalue
Fstat = res["CaseControl2","F value"]
pval = res["CaseControl2","Pr(>F)"]
etasq_res = etasq(mod2use)
aovres$df1_all[i] = res["CaseControl2","Df"]
aovres$df2_all[i] = res["Residuals","Df"]
aovres$Fstat_all[i] = Fstat
aovres$pval_all[i] = pval
aovres$etasq_all[i] = etasq_res["CaseControl2","Partial eta^2"]
# get residual for effect size
# remove variation from covariate
covname2use = c("sexM","scan_age")
beta1 = mod2use$coefficients[covname2use, drop = FALSE]
beta1[is.na(beta1)] = 0
full_model = model.matrix(~0+as.factor(CaseControl2) + as.factor(sex) + scan_age, data=D)
colnames(full_model) = c("ASD","LD_DD","TD","TypSibASD","sex","scan_age")
covname2use = c("sex","scan_age")
D$covadj = as.numeric(t(D[,y_var] - beta1 %*% t(full_model[,covname2use])))
# specific ASD pairwise comparisons
# TD vs ASD
pw_comp_res = t.test(D[D$CaseControl2=="TD",y_var],D[D$CaseControl2=="ASD",y_var])
aovres$TD_vs_ASD_t[i] = pw_comp_res$statistic
aovres$TD_vs_ASD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="TD" | D$CaseControl2=="ASD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$TD_vs_ASD_d[i] = dres$estimate
# TypSibASD vs ASD
pw_comp_res = t.test(D[D$CaseControl2=="TypSibASD",y_var],D[D$CaseControl2=="ASD",y_var])
aovres$TypSibASD_vs_ASD_t[i] = pw_comp_res$statistic
aovres$TypSibASD_vs_ASD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="TypSibASD" | D$CaseControl2=="ASD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$TypSibASD_vs_ASD_d[i] = dres$estimate
# LD_DD vs ASD
pw_comp_res = t.test(D[D$CaseControl2=="LD_DD",y_var],D[D$CaseControl2=="ASD",y_var])
aovres$LDDD_vs_ASD_t[i] = pw_comp_res$statistic
aovres$LDDD_vs_ASD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="LD_DD" | D$CaseControl2=="ASD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$LDDD_vs_ASD_d[i] = dres$estimate
# TD vs TypSibASD
pw_comp_res = t.test(D[D$CaseControl2=="TD",y_var],D[D$CaseControl2=="TypSibASD",y_var])
aovres$TD_vs_TypSibASD_t[i] = pw_comp_res$statistic
aovres$TD_vs_TypSibASD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="TD" | D$CaseControl2=="TypSibASD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$TD_vs_TypSibASD_d[i] = dres$estimate
# TD vs LD_DD
pw_comp_res = t.test(D[D$CaseControl2=="TD",y_var],D[D$CaseControl2=="LD_DD",y_var])
aovres$TD_vs_LDDD_t[i] = pw_comp_res$statistic
aovres$TD_vs_LDDD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="TD" | D$CaseControl2=="LD_DD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$TD_vs_LDDD_d[i] = dres$estimate
# TypSibASD vs LD_DD
pw_comp_res = t.test(D[D$CaseControl2=="TypSibASD",y_var],D[D$CaseControl2=="LD_DD",y_var])
aovres$TypSibASD_vs_LDDD_t[i] = pw_comp_res$statistic
aovres$TypSibASD_vs_LDDD_p[i] = pw_comp_res$p.value
Dsubset = subset(D,D$CaseControl2=="TypSibASD" | D$CaseControl2=="LD_DD")
Dsubset$CaseControl2 = factor(Dsubset$CaseControl2)
dres = effsize::cohen.d(Dsubset$covadj,Dsubset[,"CaseControl2"])
aovres$TypSibASD_vs_LDDD_d[i] = dres$estimate
# subtype model with data from TD, TD ASDSib, LD/DD, GeoASD and nonGeoASD
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"subgrp","sex","scan_age"))
subtype_model = eval(substitute(lm(formula = lm_formula, data = Dexp2_all, na.action = na.omit)))
# run ANOVA
res = anova(subtype_model)
# extract F-stat and pvalue
fstat = res["subgrp","F value"]
pval = res["subgrp","Pr(>F)"]
etasq_res = etasq(subtype_model)
aovres$df1_subtype[i] = res["subgrp","Df"]
aovres$df2_subtype[i] = res["Residuals","Df"]
aovres$fstat_subtype[i] = fstat
aovres$pval_subtype[i] = pval
aovres$etasq_subtype[i] = etasq_res["subgrp","Partial eta^2"]
# construct linear model for subtype model with meanDVARSwavelet as covariate
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s + %s",y_var,"subgrp","sex","scan_age","meanDVARSwavelet"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dtmp, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
# extract F-stat and pvalue
fstat = res[1,4]
pval = res[1,5]
aovres$fstat_subtypeDVARScov[i] = fstat
aovres$pval_subtypeDVARScov[i] = pval
# compare model with no ASD subtyping to ASD subtype model
lm_formula = as.formula(sprintf("%s ~ %s + %s +%s",y_var,"CaseControl2","sex","scan_age"))
notsubtype_model = eval(substitute(lm(formula = lm_formula, data = Dexp2_all, na.action = na.omit)))
lm_formula = as.formula(sprintf("%s ~ %s + %s +%s",y_var,"subgrp","sex","scan_age"))
subtype_model = eval(substitute(lm(formula = lm_formula, data = Dexp2_all, na.action = na.omit)))
aovres$AIC_nosubtype[i] = AIC(notsubtype_model)
aovres$AIC_subtype[i] = AIC(subtype_model)
if (aovres$AIC_subtype[i]<aovres$AIC_nosubtype[i]){
aovres$AIC_delta[i] = aovres$AIC_nosubtype[i] - aovres$AIC_subtype[i]
} else {
aovres$AIC_delta[i] = aovres$AIC_subtype[i] - aovres$AIC_nosubtype[i]
}
}
aovres$fdr_all = p.adjust(aovres$pval_all,method="fdr")
aovres$fdr_subtype = p.adjust(aovres$pval_subtype,method="fdr")
aovres$fdr_subtypeDVARScov = p.adjust(aovres$pval_subtypeDVARScov,method="fdr")
write.csv(aovres,file = file.path(resultdir,"casectrl_subtype_allcomps_output.csv"))
aovres[order(aovres$fdr_all),]## df1_all df2_all Fstat_all pval_all etasq_all fdr_all
## IC02_IC10 3 189 5.98793605 0.0006411647 0.088751689 0.01030673
## IC05_IC10 3 189 6.62772145 0.0002793470 0.099384309 0.01030673
## IC09_IC10 3 189 5.93476452 0.0006871150 0.087604923 0.01030673
## IC04_IC28 3 189 3.59599439 0.0146332362 0.052651846 0.16462391
## IC02_IC26 3 189 3.18349890 0.0250784759 0.043636902 0.17123961
## IC05_IC21 3 189 3.21542015 0.0240560685 0.049726859 0.17123961
## IC05_IC26 3 189 3.13723380 0.0266372722 0.042710392 0.17123961
## IC09_IC21 3 189 2.52563127 0.0588964755 0.045368141 0.33129267
## IC02_IC05 3 189 1.85658005 0.1384172660 0.021808414 0.41475529
## IC05_IC06 3 189 2.19003321 0.0906344498 0.030689382 0.41475529
## IC05_IC09 3 189 1.84314499 0.1407807151 0.029416041 0.41475529
## IC06_IC11 3 189 1.93596121 0.1252125277 0.015507204 0.41475529
## IC06_IC26 3 189 1.84836746 0.1398574293 0.032805613 0.41475529
## IC06_IC28 3 189 1.87143126 0.1358490509 0.024799059 0.41475529
## IC09_IC28 3 189 1.80626863 0.1474685488 0.024874876 0.41475529
## IC11_IC21 3 189 1.82247379 0.1444929199 0.027690262 0.41475529
## IC04_IC06 3 189 1.74098571 0.1600608309 0.027834059 0.42237866
## IC04_IC11 3 189 1.69778394 0.1689514649 0.020529232 0.42237866
## IC09_IC11 3 189 1.51685457 0.2115341143 0.021710412 0.50100185
## IC21_IC26 3 189 1.44706767 0.2305108843 0.021374848 0.51864949
## IC02_IC21 3 189 1.29807422 0.2764377752 0.028022073 0.59236666
## IC02_IC06 3 189 1.03640431 0.3776502418 0.014383703 0.63158570
## IC04_IC26 3 189 1.01936383 0.3852538618 0.021416531 0.63158570
## IC05_IC11 3 189 1.19323168 0.3136269284 0.014316538 0.63158570
## IC05_IC28 3 189 1.00232924 0.3929866559 0.013045065 0.63158570
## IC06_IC10 3 189 1.05942362 0.3675859303 0.011644740 0.63158570
## IC10_IC11 3 189 1.13857985 0.3347504230 0.017730870 0.63158570
## IC10_IC28 3 189 1.11609063 0.3438048973 0.015417061 0.63158570
## IC10_IC26 3 189 0.89850588 0.4430328726 0.014796062 0.68746480
## IC26_IC28 3 189 0.82825125 0.4798131092 0.016240274 0.71971966
## IC02_IC11 3 189 0.65685439 0.5795915116 0.008983540 0.73694497
## IC02_IC28 3 189 0.77557776 0.5089620687 0.015806378 0.73694497
## IC04_IC09 3 189 0.71724663 0.5428175717 0.006851961 0.73694497
## IC06_IC21 3 189 0.65039133 0.5836287384 0.006872911 0.73694497
## IC11_IC28 3 189 0.64095814 0.5895559748 0.010929312 0.73694497
## IC21_IC28 3 189 0.67193953 0.5702442027 0.012166235 0.73694497
## IC04_IC10 3 189 0.43773028 0.7262613162 0.006929300 0.85798570
## IC09_IC26 3 189 0.41337327 0.7435876074 0.004132000 0.85798570
## IC11_IC26 3 189 0.42425285 0.7358308595 0.006405970 0.85798570
## IC02_IC09 3 189 0.37978526 0.7676849435 0.007936667 0.86364556
## IC10_IC21 3 189 0.29750106 0.8271765680 0.002378506 0.90787672
## IC02_IC04 3 189 0.24917578 0.8618547779 0.001682733 0.91181365
## IC04_IC21 3 189 0.20687689 0.8915511274 0.004277508 0.91181365
## IC06_IC09 3 189 0.21549452 0.8855740320 0.001069455 0.91181365
## IC04_IC05 3 189 0.07011502 0.9758300655 0.001092861 0.97583007
## df1_subtype df2_subtype fstat_subtype pval_subtype etasq_subtype
## IC02_IC10 4 157 5.2464201 0.0005418251 0.118932212
## IC05_IC10 4 157 5.5304785 0.0003428665 0.126778502
## IC09_IC10 4 157 3.8809450 0.0049204069 0.088662416
## IC04_IC28 4 157 2.5034963 0.0445057263 0.060818378
## IC02_IC26 4 157 1.9292006 0.1081992108 0.043753148
## IC05_IC21 4 157 3.3536872 0.0115052541 0.077341312
## IC05_IC26 4 157 2.7126302 0.0320040841 0.060404262
## IC09_IC21 4 157 1.8141057 0.1287429395 0.047637578
## IC02_IC05 4 157 1.7458510 0.1426004865 0.030382788
## IC05_IC06 4 157 1.9131827 0.1108603245 0.043316172
## IC05_IC09 4 157 1.6319616 0.1688513862 0.041163129
## IC06_IC11 4 157 1.1265917 0.3459943087 0.015881961
## IC06_IC26 4 157 1.1999246 0.3130703643 0.034909283
## IC06_IC28 4 157 0.8208392 0.5136849112 0.017445395
## IC09_IC28 4 157 2.1311358 0.0794498576 0.048756673
## IC11_IC21 4 157 1.5115729 0.2013757077 0.038283935
## IC04_IC06 4 157 1.4016006 0.2359325823 0.035649341
## IC04_IC11 4 157 1.0359784 0.3905237257 0.021072844
## IC09_IC11 4 157 1.7578038 0.1400772618 0.043013462
## IC21_IC26 4 157 0.8969769 0.4673138878 0.022531226
## IC02_IC21 4 157 1.6080595 0.1748983510 0.038837899
## IC02_IC06 4 157 1.5217349 0.1984247461 0.033316960
## IC04_IC26 4 157 0.9068775 0.4615046722 0.028550142
## IC05_IC11 4 157 1.5935604 0.1786626285 0.037748710
## IC05_IC28 4 157 0.8854963 0.4741145299 0.018885306
## IC06_IC10 4 157 1.0853656 0.3657184250 0.022538450
## IC10_IC11 4 157 1.7382713 0.1442224454 0.042533110
## IC10_IC28 4 157 0.6966364 0.5953562714 0.016628342
## IC10_IC26 4 157 0.7333331 0.5705063578 0.019385136
## IC26_IC28 4 157 0.9858475 0.4170295825 0.026408372
## IC02_IC11 4 157 1.2798160 0.2802431531 0.028770090
## IC02_IC28 4 157 0.7576068 0.5543875447 0.024800645
## IC04_IC09 4 157 0.7444058 0.5631210148 0.009891132
## IC06_IC21 4 157 0.5585791 0.6930588063 0.010575431
## IC11_IC28 4 157 0.3923445 0.8139063826 0.011003304
## IC21_IC28 4 157 0.5311015 0.7130482019 0.015920087
## IC04_IC10 4 157 0.3050935 0.8742362886 0.007325977
## IC09_IC26 4 157 0.7120184 0.5848718872 0.011768388
## IC11_IC26 4 157 0.3430069 0.8485742004 0.008254619
## IC02_IC09 4 157 0.2612827 0.9023892472 0.008277410
## IC10_IC21 4 157 0.1634915 0.9565534730 0.004481284
## IC02_IC04 4 157 0.5549054 0.6957250661 0.010238600
## IC04_IC21 4 157 0.1650124 0.9558304150 0.003947705
## IC06_IC09 4 157 0.3525796 0.8419446820 0.006685714
## IC04_IC05 4 157 0.4283900 0.7879751850 0.009811804
## fdr_subtype fstat_subtypeDVARScov pval_subtypeDVARScov
## IC02_IC10 0.01219106 5.2814873 0.0005136610
## IC05_IC10 0.01219106 5.5355006 0.0003412898
## IC09_IC10 0.07380610 4.3649282 0.0022554104
## IC04_IC28 0.33379295 2.5175334 0.0435626902
## IC02_IC26 0.49923154 1.9208336 0.1096178143
## IC05_IC21 0.12943411 3.3583690 0.0114329436
## IC05_IC26 0.28803676 2.7262390 0.0313458541
## IC09_IC21 0.49923154 1.8075170 0.1300603941
## IC02_IC05 0.49923154 1.7588058 0.1399046368
## IC05_IC06 0.49923154 1.9203919 0.1096912562
## IC05_IC09 0.50248864 1.6216400 0.1714754113
## IC06_IC11 0.70771563 1.1348350 0.3421831974
## IC06_IC26 0.67086507 1.2018144 0.3122869006
## IC06_IC28 0.78797154 0.8190008 0.5148542412
## IC09_IC28 0.49923154 2.1176564 0.0811498802
## IC11_IC21 0.50343927 1.5200899 0.1989356930
## IC04_IC06 0.55878770 1.3926800 0.2389890735
## IC04_IC11 0.73223199 1.0593019 0.3786723026
## IC09_IC11 0.49923154 1.7560695 0.1404777629
## IC21_IC26 0.76196978 0.8922855 0.4701012885
## IC02_IC21 0.50248864 1.6178789 0.1724271902
## IC02_IC06 0.50343927 1.5120826 0.2012625685
## IC04_IC26 0.76196978 0.9103165 0.4595164643
## IC05_IC11 0.50248864 1.5862200 0.1806330476
## IC05_IC28 0.76196978 0.9015612 0.4646347662
## IC06_IC10 0.71553605 1.0801808 0.3682868205
## IC10_IC11 0.49923154 1.7274929 0.1465954745
## IC10_IC28 0.78797154 0.6942824 0.5969766807
## IC10_IC26 0.78797154 0.7301442 0.5726523674
## IC26_IC28 0.75065325 0.9795703 0.4204638333
## IC02_IC11 0.63054709 1.2760942 0.2817355009
## IC02_IC28 0.78797154 0.7752356 0.5428601305
## IC04_IC09 0.78797154 0.7417156 0.5649200816
## IC06_IC21 0.86722079 0.5569532 0.6942404785
## IC11_IC28 0.93136193 0.4146984 0.7978669794
## IC21_IC28 0.86722079 0.5277312 0.7155072947
## IC04_IC10 0.93668174 0.3064107 0.8733595741
## IC09_IC26 0.78797154 0.7346309 0.5696472893
## IC11_IC26 0.93136193 0.3870390 0.8176849288
## IC02_IC09 0.94436084 0.2622887 0.9017614498
## IC10_IC21 0.95655347 0.1654180 0.9556346967
## IC02_IC04 0.86722079 0.5865714 0.6728285255
## IC04_IC21 0.95655347 0.1687559 0.9540297381
## IC06_IC09 0.93136193 0.3617354 0.8355532148
## IC04_IC05 0.93136193 0.4270739 0.7889264394
## fdr_subtypeDVARScov AIC_nosubtype AIC_subtype AIC_delta
## IC02_IC10 0.01155737 -47.601616 -51.329649 3.7280333
## IC05_IC10 0.01155737 -89.543146 -87.727228 1.8159180
## IC09_IC10 0.03383116 19.207204 20.297832 1.0906282
## IC04_IC28 0.32672018 -106.458925 -104.742943 1.7159814
## IC02_IC26 0.50315642 -169.978482 -168.124064 1.8544183
## IC05_IC21 0.12862062 -78.592384 -79.526365 0.9339807
## IC05_IC26 0.28211269 -255.830825 -254.253969 1.5768561
## IC09_IC21 0.50315642 -60.658256 -58.686249 1.9720067
## IC02_IC05 0.50315642 82.308288 82.882327 0.5740396
## IC05_IC06 0.50315642 95.057050 97.054316 1.9972661
## IC05_IC09 0.50315642 -52.451869 -50.975913 1.4759564
## IC06_IC11 0.69992018 120.108604 122.103260 1.9946554
## IC06_IC26 0.66918622 -223.189912 -222.636215 0.5536971
## IC06_IC28 0.78088959 -85.724944 -83.838817 1.8861263
## IC09_IC28 0.50315642 -42.049236 -40.331381 1.7178552
## IC11_IC21 0.50315642 -57.601587 -56.778775 0.8228122
## IC04_IC06 0.56602675 -22.458550 -20.787183 1.6713667
## IC04_IC11 0.71001057 2.124502 2.899205 0.7747028
## IC09_IC11 0.50315642 72.441559 74.136379 1.6948199
## IC21_IC26 0.75551993 -254.287124 -252.355442 1.9316818
## IC02_IC21 0.50315642 39.354436 38.471250 0.8831861
## IC02_IC06 0.50315642 148.702901 147.494989 1.2079123
## IC04_IC26 0.75551993 -236.735500 -235.222713 1.5127872
## IC05_IC11 0.50315642 99.736313 99.372038 0.3642751
## IC05_IC28 0.75551993 -100.826639 -98.978566 1.8480723
## IC06_IC10 0.71001057 -84.117563 -83.221939 0.8956241
## IC10_IC11 0.50315642 -8.805833 -12.076756 3.2709235
## IC10_IC28 0.79011619 -22.210481 -20.739908 1.4705726
## IC10_IC26 0.78088959 -266.070825 -264.149009 1.9218155
## IC26_IC28 0.75551993 -217.477024 -216.549457 0.9275674
## IC02_IC11 0.63390488 106.754312 107.018529 0.2642172
## IC02_IC28 0.78088959 2.609779 4.339769 1.7299894
## IC04_IC09 0.78088959 -73.155964 -71.286644 1.8693194
## IC06_IC21 0.86780060 -65.711022 -63.886750 1.8242726
## IC11_IC28 0.91707060 4.601489 6.592426 1.9909366
## IC21_IC28 0.87021157 -52.118375 -50.233033 1.8853422
## IC04_IC10 0.93574240 -35.998229 -34.053934 1.9442944
## IC09_IC26 0.78088959 -137.025822 -136.298689 0.7271324
## IC11_IC26 0.91707060 -143.994361 -142.248185 1.7461754
## IC02_IC09 0.94370384 122.941268 124.940886 1.9996182
## IC10_IC21 0.95563470 -52.959008 -51.248201 1.7108066
## IC02_IC04 0.86506525 3.619161 4.291625 0.6724636
## IC04_IC21 0.95563470 -35.278770 -33.642334 1.6364354
## IC06_IC09 0.91707060 13.332289 14.511621 1.1793315
## IC04_IC05 0.91707060 -61.822138 -60.855373 0.9667645
## TD_vs_ASD_t TD_vs_ASD_p TD_vs_ASD_d TypSibASD_vs_ASD_t
## IC02_IC10 4.08833950 7.685064e-05 -0.661803649 2.52891112
## IC05_IC10 2.95087394 3.934468e-03 -0.521973563 2.08056846
## IC09_IC10 3.60817859 4.926817e-04 -0.662503945 1.60514818
## IC04_IC28 -1.76606802 8.025688e-02 0.309871665 -1.24243412
## IC02_IC26 0.39331789 6.949485e-01 -0.048477894 3.25085384
## IC05_IC21 2.18160767 3.191754e-02 -0.426570460 -0.93556094
## IC05_IC26 -0.04949474 9.606415e-01 0.007903278 -1.02296194
## IC09_IC21 1.12775990 2.617336e-01 -0.209555217 0.89907726
## IC02_IC05 -0.76411785 4.466252e-01 0.188520525 -1.50390486
## IC05_IC06 1.90562181 5.906742e-02 -0.295264273 -0.36031634
## IC05_IC09 -0.69732069 4.872455e-01 0.106570575 1.86430163
## IC06_IC11 -1.87524568 6.332053e-02 0.264260958 -2.12850956
## IC06_IC26 1.90066040 6.022430e-02 -0.353510307 1.75405272
## IC06_IC28 -1.81547459 7.227499e-02 0.303244712 -0.47425973
## IC09_IC28 -2.34545668 2.045466e-02 0.358652502 -0.88269113
## IC11_IC21 -0.46396030 6.437510e-01 0.074677776 0.62854280
## IC04_IC06 -1.20258706 2.315526e-01 0.219546493 -1.98821636
## IC04_IC11 1.14516866 2.548468e-01 -0.203050073 1.08230532
## IC09_IC11 1.91670998 5.815231e-02 -0.318727711 1.08027707
## IC21_IC26 -0.44160253 6.597511e-01 0.055661163 0.55691110
## IC02_IC21 1.91073845 5.852409e-02 -0.393557606 -0.07518491
## IC02_IC06 0.68804259 4.929031e-01 -0.127182434 -0.67397697
## IC04_IC26 1.25194125 2.134916e-01 -0.250988823 1.11379352
## IC05_IC11 -1.62316358 1.073264e-01 0.241994646 -1.71089498
## IC05_IC28 -1.35122981 1.795440e-01 0.199554025 -0.07819236
## IC06_IC10 0.49466501 6.217695e-01 -0.033272182 1.51861052
## IC10_IC11 1.34484249 1.815175e-01 -0.228129641 0.88711047
## IC10_IC28 0.43487622 6.644434e-01 -0.078969602 1.65381179
## IC10_IC26 1.34616600 1.813330e-01 -0.230728690 0.22047367
## IC26_IC28 1.42631969 1.571030e-01 -0.283825039 0.45501387
## IC02_IC11 -0.51021690 6.110792e-01 0.058197487 -0.74006594
## IC02_IC28 0.91101944 3.642256e-01 -0.148038199 -0.15265949
## IC04_IC09 -0.38594411 7.004998e-01 0.035614619 -1.22668061
## IC06_IC21 -0.81767650 4.155345e-01 0.143788690 -0.18806540
## IC11_IC28 1.34252287 1.821538e-01 -0.232978904 -0.14946364
## IC21_IC28 -1.30704894 1.944823e-01 0.245191838 -0.46818004
## IC04_IC10 0.60462558 5.467810e-01 -0.145824608 -0.59591614
## IC09_IC26 0.54436635 5.871097e-01 -0.051419753 0.91273213
## IC11_IC26 -0.64063517 5.231966e-01 0.115805958 0.37952573
## IC02_IC09 -0.72167774 4.721889e-01 0.153966077 -0.53188793
## IC10_IC21 -0.73670077 4.628613e-01 0.090885144 -0.73907481
## IC02_IC04 -0.27227874 7.860045e-01 -0.002581072 -0.97387047
## IC04_IC21 0.37934618 7.051324e-01 -0.088797337 0.42029412
## IC06_IC09 -0.32474611 7.458428e-01 0.060457096 0.44989804
## IC04_IC05 -0.16933350 8.658639e-01 0.045355647 0.26692174
## TypSibASD_vs_ASD_p TypSibASD_vs_ASD_d LDDD_vs_ASD_t
## IC02_IC10 0.018893459 -0.60549851 1.664981303
## IC05_IC10 0.051204745 -0.64407807 2.474492564
## IC09_IC10 0.122955211 -0.40768683 2.276207344
## IC04_IC28 0.230040294 0.40910237 -2.410397306
## IC02_IC26 0.004109426 -0.85183487 0.738127855
## IC05_IC21 0.361434072 0.27613886 1.727497912
## IC05_IC26 0.317428302 0.21264252 -2.444897337
## IC09_IC21 0.377664368 -0.28308118 2.504603484
## IC02_IC05 0.150099443 0.48781992 -1.340594150
## IC05_IC06 0.722658809 0.18702522 2.172674299
## IC05_IC09 0.077204675 -0.52389168 -1.273976165
## IC06_IC11 0.043376527 0.26219332 -1.010866622
## IC06_IC26 0.094077896 -0.47822136 0.064676883
## IC06_IC28 0.640219736 0.06354073 -1.750154604
## IC09_IC28 0.387638366 0.17989116 -1.102535802
## IC11_IC21 0.537709201 -0.22043316 1.822353970
## IC04_IC06 0.061297514 0.56384079 -0.070216971
## IC04_IC11 0.292688188 -0.27258581 1.957693543
## IC09_IC11 0.290988913 -0.19888308 0.034663185
## IC21_IC26 0.583504369 -0.19959086 -1.428512006
## IC02_IC21 0.940753947 -0.15054231 0.261587360
## IC02_IC06 0.507980049 0.19522240 1.278363544
## IC04_IC26 0.280217553 -0.44176580 -0.101391205
## IC05_IC11 0.098101467 0.27463639 -0.716777911
## IC05_IC28 0.938489913 -0.02498218 0.671441994
## IC06_IC10 0.145240410 -0.37994367 -0.334485807
## IC10_IC11 0.385761083 -0.24973969 1.674121074
## IC10_IC28 0.113873588 -0.43104971 0.889749018
## IC10_IC26 0.827123262 -0.06648209 0.933366138
## IC26_IC28 0.653949571 -0.15740136 0.428707567
## IC02_IC11 0.468298368 0.15798090 -1.402569634
## IC02_IC28 0.880346519 0.09191791 -1.141456408
## IC04_IC09 0.236284374 0.29590365 0.112568695
## IC06_IC21 0.852880588 0.13267460 1.328837538
## IC11_IC28 0.882769719 0.04765042 0.394027054
## IC21_IC28 0.644552527 0.09207453 0.012233928
## IC04_IC10 0.558359344 0.10393199 0.536532437
## IC09_IC26 0.372001463 -0.22115212 -0.295386495
## IC11_IC26 0.708822625 -0.15234046 -0.756973721
## IC02_IC09 0.601581155 0.23834872 -0.661078783
## IC10_IC21 0.468125748 0.11862669 -0.033595333
## IC02_IC04 0.341251748 0.16008024 -0.007866081
## IC04_IC21 0.679001987 -0.16532429 0.749491351
## IC06_IC09 0.656759540 -0.01121969 0.624423151
## IC04_IC05 0.792351858 -0.04112908 0.272581163
## LDDD_vs_ASD_p LDDD_vs_ASD_d TD_vs_TypSibASD_t TD_vs_TypSibASD_p
## IC02_IC10 0.11421669 -0.513481887 0.21820131 0.82899786
## IC05_IC10 0.02550088 -1.055820607 -0.25473259 0.80105022
## IC09_IC10 0.03574124 -0.633627226 1.09186292 0.28280785
## IC04_IC28 0.02814014 0.815223789 0.31394233 0.75661920
## IC02_IC26 0.47160179 -0.308303306 -2.73761457 0.01085497
## IC05_IC21 0.10102532 -0.422175267 2.17545566 0.03785921
## IC05_IC26 0.02648656 0.848796230 0.80832898 0.42360484
## IC09_IC21 0.02219477 -0.769672822 -0.05378268 0.95747040
## IC02_IC05 0.19833317 0.326487706 1.03012554 0.31402864
## IC05_IC06 0.04145219 -0.400138990 1.30414213 0.20652086
## IC05_IC09 0.21551330 0.258303636 -2.13723402 0.04174990
## IC06_IC11 0.32440801 0.106926074 0.58562082 0.56218451
## IC06_IC26 0.94913575 -0.019537513 -0.34455705 0.73288266
## IC06_IC28 0.09738586 0.431471170 -0.72520080 0.47443052
## IC09_IC28 0.28509393 0.255145793 -0.43654949 0.66680186
## IC11_IC21 0.08665476 -0.610152495 -0.82499411 0.41802720
## IC04_IC06 0.94483857 -0.043010476 1.27589466 0.21537305
## IC04_IC11 0.06616865 -0.484584063 -0.35830142 0.72318772
## IC09_IC11 0.97268223 0.019557453 0.46136210 0.64733747
## IC21_IC26 0.17273049 0.501106258 -0.79442073 0.43315542
## IC02_IC21 0.79675794 -0.151059255 1.30704722 0.20229281
## IC02_IC06 0.21826748 -0.331443623 1.04788126 0.30450639
## IC04_IC26 0.92027720 0.003682787 -0.43704057 0.66641475
## IC05_IC11 0.48249312 0.180912221 0.25743042 0.79824681
## IC05_IC28 0.51123661 -0.220245614 -0.68881003 0.49755147
## IC06_IC10 0.74157201 0.050041350 -1.19681625 0.24410639
## IC10_IC11 0.10918760 -0.392026709 -0.04835369 0.96182447
## IC10_IC28 0.38570354 -0.231220996 -1.34257728 0.19251562
## IC10_IC26 0.36394525 -0.334417336 0.91530710 0.36496532
## IC26_IC28 0.67185522 -0.058327783 0.52254239 0.60512798
## IC02_IC11 0.17703191 0.369359246 0.37745854 0.70888506
## IC02_IC28 0.26734552 0.348597433 0.59589742 0.55765203
## IC04_IC09 0.91171295 -0.119072715 0.92423307 0.36459961
## IC06_IC21 0.19537991 -0.137007642 -0.27418605 0.78629469
## IC11_IC28 0.69862137 -0.091413704 0.86424838 0.39660802
## IC21_IC28 0.99039110 -0.064479421 -0.47610974 0.63718566
## IC04_IC10 0.59871224 -0.169574791 0.89441218 0.38001646
## IC09_IC26 0.77139183 0.050474737 -0.58794193 0.56257756
## IC11_IC26 0.45911224 0.160963317 -0.67218032 0.50873015
## IC02_IC09 0.51775885 0.227523961 0.17462513 0.86307716
## IC10_IC21 0.97352880 -0.012527263 0.23486745 0.81614577
## IC02_IC04 0.99381778 0.015509251 0.67900514 0.50205163
## IC04_IC21 0.46241192 -0.184969188 -0.20261049 0.84130434
## IC06_IC09 0.53874527 -0.037002235 -0.66544372 0.51211048
## IC04_IC05 0.78831999 -0.072349239 -0.34791064 0.73084647
## TD_vs_TypSibASD_d TD_vs_LDDD_t TD_vs_LDDD_p TD_vs_LDDD_d
## IC02_IC10 0.045881596 0.27178664 0.78883297 0.09952879
## IC05_IC10 -0.111054478 -1.31878303 0.20469367 -0.40106368
## IC09_IC10 0.276270073 -0.05346514 0.95781677 0.04706909
## IC04_IC28 0.105710807 1.52486475 0.14402146 0.42206586
## IC02_IC26 -0.728931363 -0.56363663 0.58018409 -0.17938025
## IC05_IC21 0.576537989 -0.06304070 0.95015113 0.06043328
## IC05_IC26 0.158275176 2.23996165 0.03585833 0.62728392
## IC09_IC21 -0.082207199 -1.80919219 0.08509324 -0.57209348
## IC02_IC05 0.305421357 0.89779818 0.37992371 0.12077041
## IC05_IC06 0.518267195 -0.82760205 0.41582500 -0.13203292
## IC05_IC09 -0.583891361 0.57721033 0.56740562 0.17678799
## IC06_IC11 -0.001113726 -0.21215939 0.83376980 -0.16108141
## IC06_IC26 -0.109945115 1.07922666 0.29105832 0.33860106
## IC06_IC28 -0.241854456 0.64537580 0.52537973 0.12234669
## IC09_IC28 -0.194753260 -0.03805877 0.97005479 -0.07446247
## IC11_IC21 -0.256199894 -1.94384435 0.06544994 -0.56856622
## IC04_IC06 0.375046605 -0.51524959 0.61229707 -0.23110192
## IC04_IC11 -0.063536896 -1.16977868 0.25405227 -0.26577922
## IC09_IC11 0.130924733 1.30553348 0.20150639 0.38972677
## IC21_IC26 -0.239325538 1.18361533 0.25165068 0.33053421
## IC02_IC21 0.246002295 0.69657988 0.49428537 0.21781471
## IC02_IC06 0.333289844 -0.87225654 0.39331415 -0.19162534
## IC04_IC26 -0.176567701 0.90950072 0.37107779 0.28663189
## IC05_IC11 0.032269707 -0.24721816 0.80702658 -0.05705982
## IC05_IC28 -0.215253825 -1.28893756 0.21250174 -0.34860217
## IC06_IC10 -0.374887932 0.62230457 0.53972513 0.10223357
## IC10_IC11 -0.019624823 -0.62613797 0.53633179 -0.18262577
## IC10_IC28 -0.386489049 -0.64064222 0.52914984 -0.15244062
## IC10_IC26 0.172640750 -0.19441887 0.84773528 -0.08748031
## IC26_IC28 0.127646675 0.80196654 0.42698966 0.31078102
## IC02_IC11 0.087226429 0.93540737 0.35776728 0.31035875
## IC02_IC28 0.238549603 1.64824439 0.11205711 0.57484394
## IC04_IC09 0.209516600 -0.32021737 0.75155506 -0.13664995
## IC06_IC21 -0.012510941 -1.79098689 0.08027450 -0.40027637
## IC11_IC28 0.286156893 0.23243656 0.81873252 0.11263063
## IC21_IC28 -0.150328998 -0.62492373 0.53913084 -0.23156097
## IC04_IC10 0.242991514 -0.21032913 0.83552166 -0.02192483
## IC09_IC26 -0.196878466 0.51931913 0.61018418 0.08211436
## IC11_IC26 -0.248182072 0.35317004 0.72727723 0.04158209
## IC02_IC09 0.080519963 0.27920145 0.78297121 0.06309760
## IC10_IC21 0.028944918 -0.45082605 0.65585406 -0.12095542
## IC02_IC04 0.151720911 -0.13336492 0.89515493 0.01539938
## IC04_IC21 -0.083115172 -0.47120202 0.64172971 -0.10866480
## IC06_IC09 -0.091000010 -0.83399185 0.41332843 -0.12547637
## IC04_IC05 -0.083602806 -0.35136145 0.72859483 -0.11459359
## TypSibASD_vs_LDDD_t TypSibASD_vs_LDDD_p TypSibASD_vs_LDDD_d
## IC02_IC10 0.10254260 0.91915506 -0.070864696
## IC05_IC10 -1.04240230 0.30805908 0.368247201
## IC09_IC10 -0.86341158 0.39568304 0.254291759
## IC04_IC28 0.99211915 0.32950883 -0.335748988
## IC02_IC26 1.06604454 0.29776164 -0.363695815
## IC05_IC21 -1.97382473 0.05802930 0.666614220
## IC05_IC26 1.65407849 0.11224080 -0.600685885
## IC09_IC21 -1.59587579 0.12285096 0.564350100
## IC02_IC05 -0.07231240 0.94285325 0.125122105
## IC05_IC06 -1.71092620 0.09864977 0.614781381
## IC05_IC09 2.46506755 0.02035081 -0.896836915
## IC06_IC11 -0.64069240 0.52710263 0.175754238
## IC06_IC26 1.20964896 0.23643065 -0.482539092
## IC06_IC28 1.10258667 0.27971819 -0.376384806
## IC09_IC28 0.27734350 0.78356546 -0.080653820
## IC11_IC21 -0.88267335 0.38468221 0.311173663
## IC04_IC06 -1.35300216 0.18673381 0.559659060
## IC04_IC11 -0.64947667 0.52114545 0.202520485
## IC09_IC11 0.79061721 0.43567808 -0.265532296
## IC21_IC26 1.57740024 0.12874581 -0.590307665
## IC02_IC21 -0.27186917 0.78789155 0.003514785
## IC02_IC06 -1.50222101 0.14453396 0.527827593
## IC04_IC26 1.01576740 0.31893831 -0.437606250
## IC05_IC11 -0.41292011 0.68338250 0.100603629
## IC05_IC28 -0.59446892 0.55706035 0.171861894
## IC06_IC10 1.47104798 0.15225810 -0.492020033
## IC10_IC11 -0.43344264 0.66798763 0.148244438
## IC10_IC28 0.46273367 0.64710360 -0.193950166
## IC10_IC26 -0.75753909 0.45728619 0.292284124
## IC26_IC28 0.12008332 0.90532677 -0.120941141
## IC02_IC11 0.41656078 0.68008795 -0.205375668
## IC02_IC28 0.60757373 0.54852963 -0.246786678
## IC04_IC09 -1.02081219 0.31581357 0.341797730
## IC06_IC21 -0.91276694 0.37089222 0.281821449
## IC11_IC28 -0.41738406 0.67962582 0.124762408
## IC21_IC28 -0.27599135 0.78497781 0.130884775
## IC04_IC10 -0.83592461 0.41024656 0.254810216
## IC09_IC26 0.82039589 0.41959711 -0.245659414
## IC11_IC26 0.81878745 0.41970588 -0.264243925
## IC02_IC09 0.07247327 0.94272303 0.012120946
## IC10_IC21 -0.56136811 0.57886401 0.151234570
## IC02_IC04 -0.59734043 0.55559652 0.138843795
## IC04_IC21 -0.16893059 0.86705838 0.015784334
## IC06_IC09 -0.15153229 0.88061039 0.036324839
## IC04_IC05 -0.01119008 0.99114883 0.031048945Write out results for significant connections
sig_res = vars2use[aovres$fdr_all<=fdr_thresh]
sig_res## [1] "IC02_IC10" "IC05_IC10" "IC09_IC10"Dexp2 = subset(D, is.element(D$subgrp,c("GeoASD","nonGeoASD","LD_DD","TypSibASD","TD")),
select = 1:ncol(D))
Dexp2$subgrp = factor(Dexp2$subgrp)
Dexp2$Case_vs_nonASD = NA
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("LD_DD","TypSibASD"))] = "nonASD"
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("GeoASD","nonGeoASD"))] = "ASD"
Dexp2$Case_vs_nonASD = factor(Dexp2$Case_vs_nonASD)
# set up comparisons to run
comp1 = c("GeoASD","nonGeoASD","LD_DD","TypSibASD")
comp2 = list(comp = c("nonGeoASD","LD_DD","TypSibASD","TD"),
comp = c("LD_DD","TypSibASD","TD"),
comp = c("TypSibASD","TD"),
comp = c("TD"))
# set up data frame for storing multiple comparison results
colnames2use = c("compName","tstat","pvalue","fdr_q","cohensd","np.W","np.pvalue","np.fdr_q")
mcompres = data.frame(matrix(nrow = 10,ncol = length(colnames2use)))
colnames(mcompres) = colnames2use
# set up stuff to plotting effect size matrix
dmat_idx = cbind(c(1,1,1,1,2,2,2,3,3,4), c(2,3,4,5,3,4,5,4,5,5))
dMat_grpLabels = c("GeoASD","nonGeoASD","LD_DD","TDSibASD","TD")
# more stuff for plots
yLimits = list(ylim = c(-0.4,1),
ylim = c(-0.6,1),
ylim = c(-0.3,1.6))
# loop over number of significant connections
for (i in 1:length(sig_res)) {
y_var = sig_res[i]
# model using only subgrp
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"subgrp","sex","scan_age"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dexp2, na.action = na.omit)))
subtype_model = mod2use
subtype_formula = lm_formula
# anova on model using only subgrp
sigaovres = anova(subtype_model)
# model using case-control status
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"CaseControl2","sex","scan_age"))
mod2use = eval(substitute(lm(formula = lm_formula, data = Dexp2, na.action = na.omit)))
cc_model = mod2use
cc_formula = lm_formula
# compare subtype vs case-control models with AIC
subtype_vs_cc_model_compeval = rbind(AIC(subtype_model),AIC(cc_model))
rownames(subtype_vs_cc_model_compeval) = c(subtype_formula,cc_formula)
colnames(subtype_vs_cc_model_compeval) = c("AIC")
# remove sex and scan age for effect size computation
covname2use = c("sexM","scan_age")
beta1 = mod2use$coefficients[covname2use, drop = FALSE]
beta1[is.na(beta1)] = 0
full_model = model.matrix(~0+as.factor(subgrp) + as.factor(sex) + scan_age, data = Dexp2)
colnames(full_model) = c("GeoASD","nonGeoASD","LD_DD","TypSibASD","TD","sex","scan_age")
covname2use = c("sex","scan_age")
Dexp2$covadj = as.numeric(t(Dexp2[,y_var] - beta1 %*% t(full_model[,covname2use])))
# multiple comparisons on subgroup
compCount = 0
for (ic1 in 1:length(comp1)) {
for (ic2 in 1:length(comp2[ic1]$comp)) {
compCount = compCount + 1
compname1 = comp1[ic1]
compname2 = comp2[ic1]$comp[ic2]
Dcomp = subset(Dexp2, Dexp2$subgrp==compname1 | Dexp2$subgrp==compname2,
select = c("subgrp",y_var,"covadj"))
Dcomp$subgrp = factor(Dcomp$subgrp)
tres = tres = t.test(Dcomp[Dcomp$subgrp==compname1,y_var],Dcomp[Dcomp$subgrp==compname2,y_var])
dres = effsize::cohen.d(Dcomp$covadj,Dcomp[,"subgrp"])
# Added Mann-Whitney U test --------------------------
npres = wilcox.test(Dcomp[Dcomp$subgrp==compname1,y_var],Dcomp[Dcomp$subgrp==compname2,y_var])
mcompres$np.W[compCount] = npres$statistic
mcompres$np.pvalue[compCount] = npres$p.value
# ----------------------------------------------------
mcompres$compName[compCount] = sprintf("%s vs %s",compname1,compname2)
mcompres$tstat[compCount] = tres$statistic
mcompres$pvalue[compCount] = tres$p.value
mcompres$cohensd[compCount] = dres$estimate
}#for (ic2 in 1:length(comp2[ic1]$comp))
}#for (ic1 in 1:length(comp1))
# compute FDR
mcompres$fdr_q = p.adjust(mcompres$pvalue,method = "fdr")
# write.csv(mcompres, file = file.path(resultdir,sprintf("mcompres_%s.csv",y_var)))
# compute FDR on non-parametric tests ------------------
mcompres$np.fdr_q = p.adjust(mcompres$np.pvalue,method = "fdr")
write.csv(mcompres, file = file.path(resultdir,sprintf("mcompres_%s.csv",y_var)))
# ------------------------------------------------------
# RDOC model - Correlation with FixGeo across all groups
rdoc_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"FixGeo","sex","scan_age"))
rdoc_model = eval(substitute(lm(formula = rdoc_formula, data = Dexp2, na.action = na.omit)))
subgrpNOGEOFIX_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"subgrp","sex","scan_age"))
subgrpNOGEOFIX_model = eval(substitute(lm(formula = subgrpNOGEOFIX_formula, data = Dexp2, na.action = na.omit)))
# compare RDOC vs subgrp models with AIC
model_compeval = rbind(AIC(subgrpNOGEOFIX_model),AIC(rdoc_model))
rownames(model_compeval) = c(subgrpNOGEOFIX_formula,rdoc_formula)
colnames(model_compeval) = c("AIC")
# print results to screen
print(sprintf("%s: ANOVA on stratified model", y_var))
print(sigaovres)
print(etasq(subtype_model))
print(sprintf("%s: Statistics for each pairwise comparison", y_var))
print(mcompres)
print(sprintf("%s: Model comparison for subtype vs case-control models", y_var))
print(subtype_vs_cc_model_compeval)
print(sprintf("%s: Model comparison for RDOC vs subtype models", y_var))
print(model_compeval)
# make effect size matrix as heatmap figure
ngrp = length(unique(Dexp2$subgrp))
dMat = data.frame(matrix(nrow = ngrp, ncol=ngrp))
dMat[diag(x = 1,nrow=ngrp,ncol=ngrp)==1] = NA
for (ires in 1:dim(dmat_idx)[1]) {
dMat[dmat_idx[ires,1],dmat_idx[ires,2]] = abs(mcompres$cohensd[ires])
dMat[dmat_idx[ires,2],dmat_idx[ires,1]] = abs(mcompres$cohensd[ires])
}# for (ires in 1:dim(dmat_idx)[1])
rownames(dMat) = dMat_grpLabels
colnames(dMat) = dMat_grpLabels
#plot the matrix as a heatmap using the labeledHeatmap function from WGCNA
# WGCNA::sizeGrWindow(10,10)
# pdf(file = file.path(plotdir,sprintf("FC_%s_effectSize_plot.pdf",y_var)))
# par(mar = c(6, 8.5, 3, 3))
WGCNA::labeledHeatmap(Matrix = dMat,
xLabels = rownames(dMat), yLabels = colnames(dMat),
ySymbols = NULL, colorLabels = FALSE,
colors = WGCNA::blueWhiteRed(50), textMatrix = round(dMat,digits=2),
setStdMargins = FALSE, cex.text = 1.5, zlim = c(0,1),
main = paste("Effect Size (Cohen's d)"))
# dev.off()
# make scatter-boxplots
D$subgrp4plot = D$subgrp
D$subgrp4plot = factor(D$subgrp4plot,levels(D$subgrp4plot)[c(1,2,6,3:5)])
colours2use = get_ggColorHue(6)
colours2use = c(colours2use[1],colours2use[5],colours2use[6],colours2use[2:4])
xLabel = "Group"
yLabel = "Connectivity (z)"
p = ggplot(data = D, aes_string(x = "subgrp4plot", y = y_var, colour = "subgrp4plot"))
p = p + geom_jitter(size = dotSize) +
geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
guides(colour = FALSE)
p = p + scale_y_continuous(limits = yLimits[i]$ylim,
breaks = round(seq(from = yLimits[i]$ylim[1],
to = yLimits[i]$ylim[2],
by = 0.1),digits=2)) +
geom_hline(yintercept = 0, linetype = 2) +
scale_colour_manual(values = colours2use) +
xlab(xLabel) + ylab(yLabel) +
theme(axis.text.x = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
ggsave(filename = file.path(plotdir,sprintf("FC_%s_subgrp_plot.pdf",y_var)))
print(p)
# Make RDOC plot
colours2use = get_ggColorHue(6)
colours2use = c(colours2use[1],colours2use[5],colours2use[2:4])
p = ggplot(data = Dexp2, aes_string(x = "FixGeo", y = y_var)) +
geom_point(data = Dexp2,aes(colour = subgrp)) + geom_smooth(method = lm)
p = p + scale_colour_manual(values = colours2use) +
xlab("Fixation Geometric (%)") + ylab("Connectivity (z)")
ggsave(filename = file.path(plotdir,sprintf("FC_%s_RDOCgeoFix_plot.pdf",y_var)))
print(p)
}# for (i in 1:length(sig_res))## [1] "IC02_IC10: ANOVA on stratified model"
## Analysis of Variance Table
##
## Response: IC02_IC10
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp 4 0.8513 0.212831 5.2464 0.0005418 ***
## sex 1 0.0277 0.027699 0.6828 0.4098801
## scan_age 1 0.0060 0.006044 0.1490 0.7000326
## Residuals 157 6.3690 0.040567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Partial eta^2
## subgrp 0.1189322119
## sex 0.0040032808
## scan_age 0.0009480308
## Residuals NA
## [1] "IC02_IC10: Statistics for each pairwise comparison"
## compName tstat pvalue fdr_q cohensd
## 1 GeoASD vs nonGeoASD -2.6543742 1.309552e-02 0.0261910354 0.72470597
## 2 GeoASD vs LD_DD -2.7516186 1.091693e-02 0.0261910354 -0.99527153
## 3 GeoASD vs TypSibASD -3.6308198 1.041578e-03 0.0052078906 -1.29532169
## 4 GeoASD vs TD -4.6631144 8.564567e-05 0.0008564567 -1.28028563
## 5 nonGeoASD vs LD_DD -1.1328690 2.713728e-01 0.3876753747 -0.37291200
## 6 nonGeoASD vs TypSibASD -1.7213301 9.640952e-02 0.1606825407 -0.47147648
## 7 nonGeoASD vs TD -2.7045624 7.878847e-03 0.0261910354 -0.51864838
## 8 LD_DD vs TypSibASD -0.1025426 9.191551e-01 0.9191550624 -0.06855427
## 9 LD_DD vs TD -0.2717866 7.888330e-01 0.9191550624 0.09283228
## 10 TypSibASD vs TD -0.2182013 8.289979e-01 0.9191550624 0.04117492
## np.W np.pvalue np.fdr_q
## 1 291 1.138846e-02 0.0304180677
## 2 59 1.520903e-02 0.0304180677
## 3 44 1.059778e-03 0.0052988919
## 4 156 9.558253e-05 0.0009558253
## 5 379 2.714698e-01 0.3878139525
## 6 377 1.425512e-01 0.2375852823
## 7 1253 1.367676e-02 0.0304180677
## 8 122 9.534141e-01 0.9534140883
## 9 399 8.523894e-01 0.9470993045
## 10 413 7.153397e-01 0.8941745667
## [1] "IC02_IC10: Model comparison for subtype vs case-control models"
## AIC
## IC02_IC10 ~ subgrp + sex + scan_age -51.32965
## IC02_IC10 ~ CaseControl2 + sex + scan_age -47.60162
## [1] "IC02_IC10: Model comparison for RDOC vs subtype models"
## AIC
## IC02_IC10 ~ subgrp + sex + scan_age -51.32965
## IC02_IC10 ~ FixGeo + sex + scan_age -46.13133
## [1] "IC05_IC10: ANOVA on stratified model"
## Analysis of Variance Table
##
## Response: IC05_IC10
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp 4 0.7188 0.179700 5.5305 0.0003429 ***
## sex 1 0.0258 0.025780 0.7934 0.3744368
## scan_age 1 0.0024 0.002413 0.0743 0.7855749
## Residuals 157 5.1014 0.032493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Partial eta^2
## subgrp 0.126778502
## sex 0.005234935
## scan_age 0.000472840
## Residuals NA
## [1] "IC05_IC10: Statistics for each pairwise comparison"
## compName tstat pvalue fdr_q cohensd np.W
## 1 GeoASD vs nonGeoASD -0.5469707 0.58940781 0.65489757 0.1334773 490
## 2 GeoASD vs LD_DD -2.6798823 0.01384141 0.04613803 -1.0068856 54
## 3 GeoASD vs TypSibASD -2.2801524 0.02996097 0.05607706 -0.8434620 72
## 4 GeoASD vs TD -2.6438644 0.01348821 0.04613803 -0.7123278 262
## 5 nonGeoASD vs LD_DD -2.6046314 0.01893174 0.04732936 -1.0613451 224
## 6 nonGeoASD vs TypSibASD -2.2636323 0.03364624 0.05607706 -0.7090093 326
## 7 nonGeoASD vs TD -3.0678711 0.00270912 0.02709120 -0.5891286 1142
## 8 LD_DD vs TypSibASD 1.0424023 0.30805908 0.38507385 0.3688638 156
## 9 LD_DD vs TD 1.3187830 0.20469367 0.29241952 -0.4003244 550
## 10 TypSibASD vs TD 0.2547326 0.80105022 0.80105022 -0.1092939 438
## np.pvalue np.fdr_q
## 1 0.945738872 0.98353038
## 2 0.008212685 0.02737562
## 3 0.035150037 0.05992341
## 4 0.014575291 0.03643823
## 5 0.001979682 0.01063810
## 6 0.035954046 0.05992341
## 7 0.002127620 0.01063810
## 8 0.162740408 0.20342551
## 9 0.049890903 0.07127272
## 10 0.983530382 0.98353038
## [1] "IC05_IC10: Model comparison for subtype vs case-control models"
## AIC
## IC05_IC10 ~ subgrp + sex + scan_age -87.72723
## IC05_IC10 ~ CaseControl2 + sex + scan_age -89.54315
## [1] "IC05_IC10: Model comparison for RDOC vs subtype models"
## AIC
## IC05_IC10 ~ subgrp + sex + scan_age -87.72723
## IC05_IC10 ~ FixGeo + sex + scan_age -74.63179
## [1] "IC09_IC10: ANOVA on stratified model"
## Analysis of Variance Table
##
## Response: IC09_IC10
## Df Sum Sq Mean Sq F value Pr(>F)
## subgrp 4 0.9746 0.24366 3.8809 0.00492 **
## sex 1 0.1426 0.14260 2.2713 0.13380
## scan_age 1 0.4806 0.48064 7.6554 0.00634 **
## Residuals 157 9.8571 0.06278
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Partial eta^2
## subgrp 0.08866242
## sex 0.01050586
## scan_age 0.04649366
## Residuals NA
## [1] "IC09_IC10: Statistics for each pairwise comparison"
## compName tstat pvalue fdr_q cohensd
## 1 GeoASD vs nonGeoASD -1.03272101 0.312762087 0.39095261 0.27318454
## 2 GeoASD vs LD_DD -2.32828059 0.027109839 0.09036613 -0.75648006
## 3 GeoASD vs TypSibASD -1.76278020 0.088748239 0.17749648 -0.64038849
## 4 GeoASD vs TD -2.85165854 0.008475410 0.04237705 -0.81099881
## 5 nonGeoASD vs LD_DD -1.93400349 0.067070406 0.16767601 -0.51051908
## 6 nonGeoASD vs TypSibASD -1.19094727 0.243811749 0.39095261 -0.31174329
## 7 nonGeoASD vs TD -2.86843486 0.004954703 0.04237705 -0.58167336
## 8 LD_DD vs TypSibASD 0.86341158 0.395683043 0.43964783 0.23353328
## 9 LD_DD vs TD 0.05346514 0.957816772 0.95781677 0.08578674
## 10 TypSibASD vs TD -1.09186292 0.282807850 0.39095261 0.36644480
## np.W np.pvalue np.fdr_q
## 1 390 0.191725895 0.31954316
## 2 60 0.017094319 0.05698106
## 3 78 0.061458212 0.15364553
## 4 227 0.003450821 0.03450821
## 5 341 0.112188540 0.22437708
## 6 418 0.337553921 0.42194240
## 7 1227 0.009117170 0.04558585
## 8 133 0.625980888 0.69553432
## 9 410 0.977162676 0.97716268
## 10 368 0.325122217 0.42194240
## [1] "IC09_IC10: Model comparison for subtype vs case-control models"
## AIC
## IC09_IC10 ~ subgrp + sex + scan_age 20.29783
## IC09_IC10 ~ CaseControl2 + sex + scan_age 19.20720
## [1] "IC09_IC10: Model comparison for RDOC vs subtype models"
## AIC
## IC09_IC10 ~ subgrp + sex + scan_age 20.29783
## IC09_IC10 ~ FixGeo + sex + scan_age 25.42071
Cross-validation of models to find best model with lowest mean squared prediction error and mean absolute percentage error
sig_res = vars2use[aovres$fdr_all<=fdr_thresh]
Dexp2 = subset(D, is.element(D$subgrp,c("GeoASD","nonGeoASD","LD_DD","TypSibASD","TD")),
select = 1:ncol(D))
Dexp2$subgrp = factor(Dexp2$subgrp)
Dexp2$Case_vs_nonASD = NA
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("LD_DD","TypSibASD"))] = "nonASD"
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("GeoASD","nonGeoASD"))] = "ASD"
Dexp2$Case_vs_nonASD = factor(Dexp2$Case_vs_nonASD)
kfolds = 5
cols2use = c("avg_mspe_subtype","avg_mspe_casectrl")
res = data.frame(matrix(nrow = length(sig_res), ncol = length(cols2use)))
colnames(res) = cols2use
rownames(res) = sig_res
subtype_mspe_res = data.frame(matrix(nrow = kfolds, ncol=length(sig_res)))
cc_mspe_res = data.frame(matrix(nrow = kfolds, ncol=length(sig_res)))
colnames(subtype_mspe_res) = sig_res
colnames(cc_mspe_res) = sig_res
subtype_mape_res = data.frame(matrix(nrow = kfolds, ncol=length(sig_res)))
cc_mape_res = data.frame(matrix(nrow = kfolds, ncol=length(sig_res)))
colnames(subtype_mape_res) = sig_res
colnames(cc_mape_res) = sig_res
# loop over number of significant connections
for (i in 1:length(sig_res)) {
y_var = sig_res[i]
# make cross-validation indices
set.seed(999)
cvind = vfold_cv(data = Dexp2, v = 5, strata = "subgrp")
for (k in 1:kfolds){
training_mask = vector(length = dim(Dexp2)[1])
training_mask[cvind$splits[[k]]$in_id] = TRUE
test_mask = !training_mask
training_data = Dexp2[training_mask,]
test_data = Dexp2[test_mask,]
# model using only subgrp
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"subgrp","sex","scan_age"))
subtype_model = eval(substitute(lm(formula = lm_formula, data = training_data, na.action = na.omit)))
predres = predict.lm(subtype_model, newdata = test_data)
# MSPE
residual_data = test_data[,y_var] - predres
sq_error = residual_data^2
subtype_mspe_res[k,sig_res[i]] = mean(sq_error, na.rm = TRUE)
# MAPE
subtype_mape_res[k,sig_res[i]] = mean(abs((test_data[,y_var] - predres)/test_data[,y_var])*100)
# model using case-control
lm_formula = as.formula(sprintf("%s ~ %s + %s + %s",y_var,"CaseControl2","sex","scan_age"))
cc_model = eval(substitute(lm(formula = lm_formula, data = training_data, na.action = na.omit)))
predres = predict.lm(cc_model, newdata = test_data)
# MSPE
residual_data = test_data[,y_var] - predres
sq_error = residual_data^2
cc_mspe_res[k,sig_res[i]] = mean(sq_error, na.rm = TRUE)
# MAPE
cc_mape_res[k,sig_res[i]] = mean(abs((test_data[,y_var] - predres)/test_data[,y_var])*100)
}#for (k in 1:kfolds){
}#for (i in 1:length(sig_res))
final_res = data.frame(rbind(colMeans(subtype_mspe_res),colMeans(cc_mspe_res)))
rownames(final_res) = c("subtype_model","casectrl_model")
write.csv(final_res,file = file.path(resultdir,"cv_mse_model_comparison.csv"))
final_res## IC02_IC10 IC05_IC10 IC09_IC10
## subtype_model 0.04194458 0.03285194 0.0632078
## casectrl_model 0.04301352 0.03264858 0.0634412final_res = data.frame(rbind(colMeans(subtype_mape_res),colMeans(cc_mape_res)))
rownames(final_res) = c("subtype_model","casectrl_model")
write.csv(final_res,file = file.path(resultdir,"cv_mape_model_comparison.csv"))
final_res## IC02_IC10 IC05_IC10 IC09_IC10
## subtype_model 125.5452 152.8099 156.0619
## casectrl_model 135.2241 152.1738 152.8081Permutation tests on specific significant connections and looking for subtype differences
nperm = 10000
set.seed(1)
sig_res = vars2use[aovres$fdr_all<=fdr_thresh]
sig_res## [1] "IC02_IC10" "IC05_IC10" "IC09_IC10"Dexp2 = subset(D, is.element(D$subgrp,c("GeoASD","nonGeoASD","LD_DD","TypSibASD","TD")),
select = 1:ncol(D))
Dexp2$subgrp = factor(Dexp2$subgrp)
Dexp2$Case_vs_nonASD = NA
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("LD_DD","TypSibASD"))] = "nonASD"
Dexp2$Case_vs_nonASD[is.element(Dexp2$subgrp,c("GeoASD","nonGeoASD"))] = "ASD"
Dexp2$Case_vs_nonASD = factor(Dexp2$Case_vs_nonASD)
# set up comparisons to run
comp1 = c("GeoASD","nonGeoASD","LD_DD","TypSibASD")
comp2 = list(comp = c("nonGeoASD","LD_DD","TypSibASD","TD"),
comp = c("LD_DD","TypSibASD","TD"),
comp = c("TypSibASD","TD"),
comp = c("TD"))
# set up data frame for storing multiple comparison results
colnames2use = c("compName","tstat")
mcompres_perm = list()
# loop over number of significant connections
for (i in 1:length(sig_res)) {
mcompres_perm[[i]] = data.frame(matrix(nrow = 10,ncol = length(colnames2use)))
colnames(mcompres_perm[[i]]) = colnames2use
mcompres_perm[[i]]$tstat = 1
for (iperm in 1:nperm){
# permute group label
Dexp2$subgrp_perm = sample(Dexp2$subgrp)
y_var = sig_res[i]
# multiple comparisons on subgroup
compCount = 0
for (ic1 in 1:length(comp1)) {
for (ic2 in 1:length(comp2[ic1]$comp)) {
compCount = compCount + 1
compname1 = comp1[ic1]
compname2 = comp2[ic1]$comp[ic2]
# statistic on real data
Dcomp_real = subset(Dexp2, Dexp2$subgrp==compname1 | Dexp2$subgrp==compname2,
select = c("subgrp",y_var))
Dcomp_real$subgrp = factor(Dcomp_real$subgrp)
tres_real = t.test(Dcomp_real[Dcomp_real$subgrp==compname1,y_var],
Dcomp_real[Dcomp_real$subgrp==compname2,y_var])
# statistic on permuted data
Dcomp_perm = subset(Dexp2, Dexp2$subgrp_perm==compname1 | Dexp2$subgrp_perm==compname2,
select = c("subgrp_perm",y_var))
Dcomp_perm$subgrp = factor(Dcomp_perm$subgrp)
tres_perm = t.test(Dcomp_perm[Dcomp_perm$subgrp_perm==compname1,y_var],
Dcomp_perm[Dcomp_perm$subgrp_perm==compname2,y_var])
# fill in mcompres_perm
mcompres_perm[[i]]$compName[compCount] = sprintf("%s vs %s",compname1,compname2)
if (abs(tres_perm$statistic) >= abs(tres_real$statistic)){
mcompres_perm[[i]]$tstat[compCount] = mcompres_perm[[i]]$tstat[compCount]+1
} # if
}#for (ic2 in 1:length(comp2[ic1]$comp))
}#for (ic1 in 1:length(comp1))
} # for (iperm in 1:nperm)
mcompres_perm[[i]]$pval = mcompres_perm[[i]]$tstat/(nperm+1)
mcompres_perm[[i]]$fdr = p.adjust(mcompres_perm[[i]]$pval, method = "fdr")
}# for (i in 1:length(sig_res))
sig_res[1]## [1] "IC02_IC10"mcompres_perm[[1]]## compName tstat pval fdr
## 1 GeoASD vs nonGeoASD 156 0.01559844 0.03119688
## 2 GeoASD vs LD_DD 106 0.01059894 0.02649735
## 3 GeoASD vs TypSibASD 17 0.00169983 0.00849915
## 4 GeoASD vs TD 4 0.00039996 0.00399960
## 5 nonGeoASD vs LD_DD 2665 0.26647335 0.38067622
## 6 nonGeoASD vs TypSibASD 939 0.09389061 0.15648435
## 7 nonGeoASD vs TD 82 0.00819918 0.02649735
## 8 LD_DD vs TypSibASD 9122 0.91210879 0.91210879
## 9 LD_DD vs TD 7878 0.78772123 0.91210879
## 10 TypSibASD vs TD 8239 0.82381762 0.91210879sig_res[2]## [1] "IC05_IC10"mcompres_perm[[2]]## compName tstat pval fdr
## 1 GeoASD vs nonGeoASD 5961 0.59604040 0.66226711
## 2 GeoASD vs LD_DD 118 0.01179882 0.03449655
## 3 GeoASD vs TypSibASD 284 0.02839716 0.04866180
## 4 GeoASD vs TD 121 0.01209879 0.03449655
## 5 nonGeoASD vs LD_DD 138 0.01379862 0.03449655
## 6 nonGeoASD vs TypSibASD 292 0.02919708 0.04866180
## 7 nonGeoASD vs TD 25 0.00249975 0.02499750
## 8 LD_DD vs TypSibASD 3102 0.31016898 0.38771123
## 9 LD_DD vs TD 2002 0.20017998 0.28597140
## 10 TypSibASD vs TD 8060 0.80591941 0.80591941sig_res[3]## [1] "IC09_IC10"mcompres_perm[[3]]## compName tstat pval fdr
## 1 GeoASD vs nonGeoASD 3148 0.31476852 0.3934607
## 2 GeoASD vs LD_DD 282 0.02819718 0.0939906
## 3 GeoASD vs TypSibASD 868 0.08679132 0.1735826
## 4 GeoASD vs TD 90 0.00899910 0.0449955
## 5 nonGeoASD vs LD_DD 631 0.06309369 0.1577342
## 6 nonGeoASD vs TypSibASD 2404 0.24037596 0.3934607
## 7 nonGeoASD vs TD 55 0.00549945 0.0449955
## 8 LD_DD vs TypSibASD 4016 0.40155984 0.4461776
## 9 LD_DD vs TD 9569 0.95680432 0.9568043
## 10 TypSibASD vs TD 2828 0.28277172 0.3934607Plot connectivity ADOS social affect relationships
IC02-IC10 ADOS Social Affect
yLimits = c(0,25)
comp2plot = "IC02_IC10"
geo_res = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
nongeo_res = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
corr_res = rbind(geo_res,nongeo_res)
rownames(corr_res) = c("GeoASD","nonGeoASD")
corr_res## r p ci95lo ci95hi
## GeoASD -0.78979061 0.00196584 -0.9958642 -0.3516584
## nonGeoASD 0.06866816 0.64455412 -0.2024230 0.3311468n = sum(Dexp2$subgrp=="GeoASD")
n2 = sum(Dexp2$subgrp=="nonGeoASD")
r_comp_res = paired.r(geo_res$r,nongeo_res$r,n=n,n2=n2)
comp_res = data.frame(matrix(nrow = 1, ncol = 2))
colnames(comp_res) = c("z","p")
rownames(comp_res) = sprintf("GeoASD vs nonGeoASD: %s, ADOS_CoSoTot",comp2plot)
comp_res$z = r_comp_res$z
comp_res$p = r_comp_res$p
comp_res## z p
## GeoASD vs nonGeoASD: IC02_IC10, ADOS_CoSoTot 3.719662 0.0001994899geo_adj = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
nongeo_adj = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
Dtmp = rbind(geo_adj,nongeo_adj)
Dconn = merge(Dfmri, Dtmp, by = "subjectId")
cols2use = get_ggColorHue(6)
cols2use = c(cols2use[1],cols2use[5])
g = ggplot(data = Dconn, aes(y = y_adj, x = x_adj, colour=subgrp2, fill = subgrp2))
g = g + geom_smooth(method = rlm, method.args = list(psi = psi.bisquare)) + geom_point() +
scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use)
yLabel = "ADOS Social Affect"
xLabel = "Connectivity Strength (z)"
g = g + xlab(xLabel) + ylab(yLabel) + theme(axis.text.x = element_text(size=fontSize+10,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
g = g + scale_y_continuous(limits = yLimits)
ggsave(filename = file.path(plotdir,sprintf("FC_%s_ADOSsoceng_corr_plot.pdf",comp2plot)))
g
IC05-IC10 ADOS Social Affect
comp2plot = "IC05_IC10"
geo_res = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
nongeo_res = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
corr_res = rbind(geo_res,nongeo_res)
rownames(corr_res) = c("GeoASD","nonGeoASD")
corr_res## r p ci95lo ci95hi
## GeoASD 0.04567706 0.8810737 -0.8157176 0.7518075
## nonGeoASD 0.22399054 0.1260338 -0.1023718 0.4920886n = sum(Dexp2$subgrp=="GeoASD")
n2 = sum(Dexp2$subgrp=="nonGeoASD")
r_comp_res = paired.r(geo_res$r,nongeo_res$r,n=n,n2=n2)
comp_res = data.frame(matrix(nrow = 1, ncol = 2))
colnames(comp_res) = c("z","p")
rownames(comp_res) = sprintf("GeoASD vs nonGeoASD: %s, ADOS_CoSoTot",comp2plot)
comp_res$z = r_comp_res$z
comp_res$p = r_comp_res$p
comp_res## z p
## GeoASD vs nonGeoASD: IC05_IC10, ADOS_CoSoTot 0.5944945 0.5521814geo_adj = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
nongeo_adj = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
Dtmp = rbind(geo_adj,nongeo_adj)
Dconn = merge(Dfmri, Dtmp, by = "subjectId")
cols2use = get_ggColorHue(6)
cols2use = c(cols2use[1],cols2use[5])
g = ggplot(data = Dconn, aes(y = y_adj, x = x_adj, colour=subgrp2, fill = subgrp2))
g = g + geom_smooth(method = rlm, method.args = list(psi = psi.bisquare)) + geom_point() +
scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use)
yLabel = "ADOS Social Affect"
xLabel = "Connectivity Strength (z)"
g = g + xlab(xLabel) + ylab(yLabel) + theme(axis.text.x = element_text(size=fontSize+10,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
g = g + scale_y_continuous(limits = yLimits)
ggsave(filename = file.path(plotdir,sprintf("FC_%s_ADOSsoceng_corr_plot.pdf",comp2plot)))
g
IC09-IC10 ADOS Social Affect
comp2plot = "IC09_IC10"
geo_res = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
nongeo_res = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_corrRes.txt",comp2plot)))
corr_res = rbind(geo_res,nongeo_res)
rownames(corr_res) = c("GeoASD","nonGeoASD")
corr_res## r p ci95lo ci95hi
## GeoASD -0.3389316 0.2563662 -0.83704174 0.5817597
## nonGeoASD 0.2950895 0.0452685 0.05118179 0.4937473n = sum(Dexp2$subgrp=="GeoASD")
n2 = sum(Dexp2$subgrp=="nonGeoASD")
r_comp_res = paired.r(geo_res$r,nongeo_res$r,n=n,n2=n2)
comp_res = data.frame(matrix(nrow = 1, ncol = 2))
colnames(comp_res) = c("z","p")
rownames(comp_res) = sprintf("GeoASD vs nonGeoASD: %s, ADOS_CoSoTot",comp2plot)
comp_res$z = r_comp_res$z
comp_res$p = r_comp_res$p
comp_res## z p
## GeoASD vs nonGeoASD: IC09_IC10, ADOS_CoSoTot 2.144411 0.03199995geo_adj = read.delim(file.path(tidydatadir,
sprintf("GeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
nongeo_adj = read.delim(file.path(tidydatadir,
sprintf("nonGeoASD_%s_ADOS_CoSoTot_xy_adjdata.txt",comp2plot)))
Dtmp = rbind(geo_adj,nongeo_adj)
Dconn = merge(Dfmri, Dtmp, by = "subjectId")
cols2use = get_ggColorHue(6)
cols2use = c(cols2use[1],cols2use[5])
g = ggplot(data = Dconn, aes(y = y_adj, x = x_adj, colour=subgrp2, fill = subgrp2))
g = g + geom_smooth(method = rlm, method.args = list(psi = psi.bisquare)) + geom_point() +
scale_colour_manual(values = cols2use) + scale_fill_manual(values = cols2use)
yLabel = "ADOS Social Affect"
xLabel = "Connectivity Strength (z)"
g = g + xlab(xLabel) + ylab(yLabel) + theme(axis.text.x = element_text(size=fontSize+10,hjust=0.5,vjust=0.5,face="plain"),
axis.text.y = element_text(size=fontSize+5,hjust=1,vjust=0,face="plain"),
axis.title.x = element_text(size=fontSize+5,hjust=0.5,vjust=0,face="plain"),
axis.title.y = element_text(size=fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
strip.text.x = element_text(size = fontSize+5,hjust=0.5,vjust=0.5,face="plain"),
plot.title = element_text(size=fontSize,hjust=0.5,vjust=0.5,face="plain"))
g = g + scale_y_continuous(limits = yLimits)
ggsave(filename = file.path(plotdir,sprintf("FC_%s_ADOSsoceng_corr_plot.pdf",comp2plot)))
g